Emergence of Chaotic Itinerancy in Simple Ecological Systems 
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Chaotic itinerancy is a universal dynamical concept that describes itinerant motion among many 
different ordered states through chaotic transition in dynamical systems. Unlike the expectation of 
the prevalence of chaotic itinerancy in high-dimensional systems, we identify chaotic itinerant behav- 
ior from a relatively simple ecological system, which consists only of two coupled consumer-resource 
pairs. The system exhibits chaotic bursting activity, in which the explosion and the shrinkage of 
the population alternate indefinitely, while the explosion of one pair co-occurs with the shrinkage of 
the other pair. We analyze successfully the bursting activity in the framework of chaotic itinerancy, 
and find that large duration times of bursts tend to cluster in time, allowing the effective burst 
prognosis. We also investigate the control schemes on the bursting activity, and demonstrate that 
invoking the competitive rise of the consumer in one pair can even elongate the burst of the other 
pair rather than shorten it. 



Since attractors determine the long-term behavior of 
dynamical systems, the concept of attractors is central 
to the analysis of many natural and artificial systems 
[H- In general, the phase space of a nonlinear dynami- 
cal system is partitioned into various basins of attraction 
from which states evolve towards the respective attrac- 
tors. These stable attractors can lose their stability with 
a change of the system condition such that the basin of 
attraction of each attractor becomes connected to each 
other through unstable manifolds. Hence, a dynamical 
state which sequentially traces out all of the destabilized 
attractor ruins emerges. This is referred to as a chaotic 
itinerant state Q. The notion of chaotic itinerancy has 
received considerable attention in studying the adapt- 
ability of complex systems, especially in relation to brain 
information processing 0, Q ■ 



To embed a chaotic itinerant state, a system is ex- 
pected to have a high degree of complexity; there- 
fore, models of chaotic itinerancy are mostly built on 
high-dimensional phase space 0]. Albeit relatively low- 
dimensional systems, two coupled Morris-Lecar neural 
oscillators were found to exhibit chaotic itinerancy Q, 
the result seems to be limited to a rather special case 
obtained by using sophisticated forms of model equa- 
tions in neurobiological systems. In the present work, 
we report that low-dimensional chaotic itinerancy ex- 
ists and arises naturally in simple ecological systems, of 
which consumer-resource dynamics has broad relevance 
in metabolic, immune, social, and economical systems. 
The wide variety of related disciplines aside, the math- 
ematical simplicity of our low-dimensional system ren- 
ders the global organization of a chaotic itinerant state 
tractable with a detailed illustration. 



At the outset, we suggest the equations of two 



consumer-resource pairs coupled via resource sharing 



dt 

dR\(2) 
dt 



aCi( 2 ) 

#1(2) - 



#1(2 



DR 



2(1) 



#1(2) 



#1(2) « 
[Cl(2) " 



R 



2(1) 



kCl(2) , 



DO 



2(1)J#1(2) 



#1(2 



(1) 



where Ci(2) and #1(2) represent the population of con- 
sumer 1 (2) and resource 1 (2), respectively, a and b de- 
note the growth and decay rates of the population of con- 
sumers, k concerns the satiability level of the consumers 
taking resources. For simplicity in our analysis, we do 
not distinguish the parameter sets of the two consumer- 
resource pairs. In this equation, #1(2) is taken by C±(2) 
primarily, as well as by C2(i) with a relative small uptake 
rate D, which ranges from to 1. 

When D is equal to zero, Eq. ([1]) splits into two Holling 
type-II forms of Lotka-Volterra equations 0, and the 
populations of each consumer-resource pair can exhibit a 
limit cycle oscillation. If D takes a nonzero value close 
to or 1, synchronous limit cycle oscillation between the 
two consumer-resource pairs arises. Complicated dynam- 
ics develop at intermediate range of D, where we can ob- 
serve irregular bursting activities as in Figs. [3a)-QJc). 
Time trajectories of C\ and C2 in Fig. QJc) show the 
bursting behaviors, and those of R\ and R2 show the 
similar patterns as well. C\ and C2, or R± and R2, fire 
bursts in an antiphase-synchronized way, such that the 
explosion of C\ or R\ co-occurs with the shrinkage of 
C2 or i?2, and vice versa. It is worth noting that other 
equations with similar systems to Eq. ([I]) also support the 
existence of such antiphase-synchronized irregular bursts 
8\. For numerical simulations, we use the parameters 
a = 2, b = 0.82, k = 0.33, and D = 0.57 unless specified. 

To check whether the apparent irregularity of the 
bursts implies their initial condition sensitiveness, we 
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FIG. 1: (a) and (b): parameter regime of irregular bursting 
activities in Eq. {TJ with a = 2, when (a) b = 0.82 or (b) k = 
0.33. (c) Time series of Ci and Ci in bursting activity. The 
arrowed time interval is magnified in (e). (d) Half-life time th 
of burst reproducibility with initial condition difference e. (e) 
Magnification of the arrowed time interval in (c). Solid lines 
denote Ci and Ri, and dotted lines C2 and Ri- 



evaluate 

M m (t) = ±J* H[C m (t')]H[C' m (t')]dt' , (2) 

where H(X) = 1 if X is in the burst mode, other- 
wise H{X) = —1. The antiphase synchronization of 
the bursts enables us to define the burst mode un- 
ambiguously, such that once Ci(2) exceeds CWi), Ci(2) 
enters the burst mode, and C*2(i) enters the shrink- 
age mode. Ci(2)(^) is calculated in the same way as 
Ci(2){t), but is initially perturbed from Ci^){t) with 
e = |q (2) (0)-C 1(2) (0)|/C 1(2) (0) « 1. Therefore, Eq. © 
gives the similarity between the bursting times of Ci(2) (t) 
and those of C[^(t) with slightly different initial condi- 
tions. If the bursting times of Ci(2)(i) and C'i( 2 )(t) are in 
complete agreement, Mu 2 \(t) = 1, whereas with no cor- 
relation between them, Mi( 2 )(t) = 0. In the following, 



we drop the subscript of Mx(2)(i) due to the statistical 
equivalence of C\{t) and C2(t). One can employ M(t) 
for determining the necessary time for the discrepancy 
to be significant between the bursting times of Ci( 2 )(t) 
and of C[( 2 j(t). It is observed that M(t) evolves rapidly 
from 1 to in the irregular bursting regime; thus, the 
half- life time of M(t) can serve as the characteristic 
time scale of the discrepancy growth. Figure [3Jd) shows 
that Th scales logarithmically to e, and using i?i(2)(i) 
and R'^Jt) instead of C 1 [ 2 ){t) and C[^(t) in Eq. ([2|) 
does not alter the current result. This logarithmic scal- 
ing reveals that the bursting is sensitive to the initial 
conditions, i.e., behaves chaotically 

To address such antiphase-synchronized chaotic bursts 
in detail, we divide a period of bursts into four stages - 
SI, SI', S2, and S2', as in Fig. [He). In stage SI, C x and 
Ri dominate C 2 and R 2 , while C\, which is supported 
primarily by R\, depresses severely the growth of R 2l 
and thereby of C 2 . Nonetheless, R 2 increases gradually 
in the negligible presence of C 2 , and R\ comes to decline 
with overpopulated C\ which takes both R\ and R 2 . In 
stage SI', the resultant shrinkage of R\ ensures that C\ 
depends mostly on R 2 for survival. Meanwhile, R 2 can 
boost the increase of C 2 , which then suppresses both R 2 
and thereby leading to the drastic decay of C\ in 
stage S2 The dominance of C 2 and R 2 in stage 

S2 is totally symmetric to that of C\ and R\ in stage SI. 
Accordingly, stage S2' analogous to stage SI' follows, and 
leads to stage SI for the next period. 

Each stage occupies finite time-span, forming a quasi- 
stable dynamical state. The alternating dominance of 
each species along the stages may be equivalent to 
the switching events among the sets of attractor ru- 
ins. In order to elucidate the underlying attractor ruin 
for a given stage, we consider an invariant subspace of 
[C\i Ri, C 2 , R 2 ) which contains only the species govern- 
ing the stage [see Fig. [Tf e)] : at stage SI, [C\, R\, 0, 0); 
at stage SI', (Ci, 0, 0, i? 2 ); at stage S2, (0, 0, C 2 , R 2 ); at 
stage S2', (0, C 2 , 0). The populations confined within 
each invariant subspace approach their own asymptotic 
solution. For instance, the limit cycle oscillation of C\ 
and R\ characterizes the asymptotic solution in the in- 
variant subspace of stage SI and thus underlies the burst- 
ing activity at this stage. Figure Ufa) shows an ac- 
tual time trajectory of the populations in phase space, 
which also embeds the asymptotic solutions in the in- 
variant subspaces. Near an invariant asymptotic solu- 
tion, the trajectory remains there for a long time, but 
finally escapes towards another invariant solution at the 
next stage. 

This escaping event is due to an existence of unstable 
manifolds outward from an invariant subspace. Along the 
transverse direction of the invariant subspace, we then 
perform a linear stability analysis to find the unstable 
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FIG. 2: (a) Projection of time trajectory onto (Ri — R2, C\ — 
Ci) for At = 3000, D = 0.57 (thin line, red in color). Over- 
lapped is invariant solution at each stage (black), (b) De- 
picted by the similar way with (a), when D — 0.4. (c) Distri- 
bution of burst duration T of Ci with stages SI and Si', or of 
Ci with stages S2 and S2', for different values of D (Ci and 
Ci show the identical distribution). Similar results are also 
observed with each of the stages. 
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SI : 
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dSR 2 
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where C\ of stage SI is evaluated in the absence of C 2 
and R 2 ■ Stages S2 and S2' take the formula obtained sim- 
ply by interchanging C\ and C 2l R\ and R 2 of stages SI 
and SI' in Eq. ([3]). The formula of unstable manifolds re- 
veals which species causes the instability of a given stage; 
the instability of stages SI and SI' is invoked by the in- 
crease of R 2 and C2, as described above in the ecologi- 
cal argument. From Eq. (|3|), we recognize that increase 
of D tends to reduce the coefficients in the right-hand 
sides, i.e., to decrease the escape rates along the unsta- 
ble manifolds. The resultant relaxation of the switching 



T 

1 n+1 



55 ~ 
45- 

35- 

- 

25- 
15 









\ 


'V 




L 

1 


\/\ 






Mr , 




>" -r 






(b) 



1 n+1 



■re 



(d) n 


35 50 65 80 


10° 






V — B 


10- 2 


\ -— ■ B' 


P(N) 




10" 4 




10" 6 





10 20 30 
N 



FIG. 3: (a) and (b): return time map for burst initiations 
when (a) D = 0.5 or (b) D = 0.57. (c) and (d): distribution 
of residual times N within boxed area in (a) [(c)] or in (b) 
[(d)]- 



events among attractor ruins is identified by comparing 
Figs. Ufa) and(2ljb), where consistently the trajectory be- 
tween the invariant solutions looks less intermingled with 
increased D. The most evident effect of elongated res- 
idence in attractor ruins appears in Fig. (He): the dis- 
tribution of burst duration T shifts to large value as D 
increases. Moreover, the larger D is, the higher the right 
peak of the bimodal duration distribution is, relatively to 
the left peak. We conclude that in the chaotic bursting 
regime, the enhanced coupling strength induces either of 
the consumer-resource pairs to dominate the other for a 
long time by an elongated bursting activity. 

To investigate the bimodality appearing in Fig. [DJc), 
we plot a return time map for burst initiations of con- 
sumers by considering the terms between the initiation of 
stage SI and that of stage S2, and between the initiation 
of stage S2 and that of stage SI, and so on. Figures[3la) 
and OJb) display nearly one-dimensional curves of such 
return time maps, where stepwise jumps between the up- 
per and lower extremes are responsible for the bimodality 
in Fig.[2fc). The relative expansion of upper extremes in 
Fig. [3jb) induces the right side of the duration distribu- 
tion in Fig. [^c) to be highly peaked. 

Referring to the return time maps, we find that map- 
ping trajectories are frequently trapped in boxed area 
B in Figs. [3ja) and[3jb). The distribution of residual 
times N (total number of iteration) within boxed area B 
is indeed more right-skewed than those of any other ar- 
eas (e.g., B') with the same size [Figs. EJc) and^d)]. 
The distribution within area B follows exponential fit 
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FIG. 4: (a) and (b): unperturbed time-series of populations. 
(c)-(f): initially the same as (a) and (b), but subjected to 
the perturbation during the arrowed time interval by reducing 
R2 — i?| m h a rf [(c), (d)] or doubling aC2 [(e), (f)], given the 
terms in Eq. {T}. 



P(N) oc p N with p = 0.64 for D = 0.5 and p = 0.68 
for D = 0.57, and shows significantly larger p than the 
surrogated data (p = 0.33 for D = 0.5, p — 0.28 for 
D = 0.57) with the same distribution of burst duration 
time. This indicates the residual dynamics within area B 
follows a poissonian process, but with considerable sur- 
vivability p per iteration. Since area B corresponds to 
relatively long durations of bursts, large duration times 
of bursts tend to cluster in time, as partially observed in 
Fig. [ljc). In this regard, the known information about 
the past duration could be beneficial to improve the burst 
prognosis efficiently. 

An important outcome of the stability analysis on at- 
tractor ruins is the application to a control scheme on 
burst duration. Since a rise of i?2(i) destabilizes the dom- 
inance of Ci(2) and Ri(2) , manually repressing the growth 
of i?2(i) might elongate the bursts of Ci(2) and i?i(2)- As 
shown in Figs.HJc) and@Jd), the repressed growth of R2 
prevents the overpopulation of C\ for a while and thus 
delays the shrinkage of R\ as well as of C\ . It should be 
noticed that the effect of delayed rise of C2 herein is not 
so essential to elongating the burst duration, despite the 
resource competition between C2 and C\. Counterintu- 
itively, even promoting the growth of C2 could be helpful 
to the burst duration, if resulting in the depression of Ri- 
Figuresllje) and[3Jf) indeed illustrate the possibility that 
a sufficiently large perturbation to increase C2 withdraws 



transiently R\ and C\, but also delays t he g rowth of R2 
thereby elongating the bursting activity 

In summary, we investigated a simple dynamical sys- 
tem, which consists only of two consumer-resource pairs 
but exhibits chaotic itinerancy naturally. The mathe- 
matical simplicity of the system gives rise to a clear view 
of the organization of a chaotic itinerant state, where 
each consumer-resource relationship underlies its corre- 
sponding attractor ruin as a dynamical 'building block'. 
Such concept of building blocks could be generically uti- 
lized when one designs other systems exhibiting chaotic 
itinerancy. In addition, analysis on a chaotic itinerant 
state was found to be applicable to the prognosis and 
control of the dynamical system, in the rather counterin- 
tuitive way. Beyond the suggested ecological system, any 
dynamical system which shows antiphase-synchronized 
chaotic bursts might be analyzable in the framework of 
chaotic itinerancy via our methodology. We expect that 
host-parasitoid systems with whiteffies and their para- 
sitic wasps could be employed for experimental valida- 
tion of our results, since parasitic wasps (consumers) are 
known to have overlapped hosts (resources) in a manner 
similar to the present model |l2j. 
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